# Title : Replication code for "Does Repression Undermine Opposition Demands? 
#         The Case of the Hong Kong National Security Law"
# Author: Tetsuro Kobayashi, Jaehyun Song, and Polly Chan
# Date  : Sep. 9, 2021

# -----------------------------------------------------------------------------#

#############################
# Load Packages and Dataset #
#############################

pacman::p_load(tidyverse, cregg, kableExtra, summarytools)

df_w1 <- read_csv("JJPS_KSC_Wave1.csv") # Wave 1
df_w2 <- read_csv("JJPS_KSC_Wave2.csv") # Wave 2

df_w1 <- df_w1 %>%
    mutate(across(A:PID, factor))

df_w2 <- df_w2 %>%
    mutate(across(A:PID, factor))

# Descriptive Statistics
# All the attribution and level names are abbreviated because some level names 
# are overlapped. See the next section to check full names.

df_w1 %>% # Wave 1
    dfSummary(style = "grid", plain.ascii = FALSE, graph.magnif = 0.85) %>%
    view()

df_w1 %>% # Wave 2
    dfSummary(style = "grid", plain.ascii = FALSE, graph.magnif = 0.85) %>%
    view()

###############
# Preparation #
###############

# Define y-axis labels
Scale_Labels = c(
    "A"  = "A commission of inquiry into alleged police brutality",
    "A1" = "Denied",
    "A2" = "Setting up an independent review committee with members appointed by the Hong Kong government",
    "A3" = "Setting up an independent commission of inquiry",
    "B"  = "Retracting the classification of protesters as rioters",
    "B1" = "Denied",
    "B2" = "Fully accepted",
    "C"  = "Amnesty for arrested protesters",
    "C1" = "Denied",
    "C2" = "Amnesty only for protestors arrested over minor offenses",
    "C3" = "Amnesty only for juvenile protestors",
    "C4" = "Fully accepted",
    "D"  = "Dual universal suffrage",
    "D1" = "Status quo",
    "D2" = "Voters may elect the CE from the candidates nominated by the Nominating Committee through 'one person, one vote'",
    "D3" = "Voters may elect all members of the Legislative Council through 'one person, one vote'",
    "D4" = "Voters may elect the CE and all members of the Legislative Council through 'one person, one vote'",
    "E"  = "Resignation of the Chief Executive",
    "E1" = "Denied",
    "E2" = "Fully accepted",
    "F"  = "Cash handout",
    "F1" = "HKD 10,000",
    "F2" = "HKD 25,000",
    "G"  = "Salary and profit tax reduction",
    "G1" = "75% of the first $20,000",
    "G2" = "100% of the first $15,000",
    "G3" = "75% of the first $40,000",
    "G4" = "100% of the first $30,000",
    "H"  = "Public housing rental exemption",
    "H1" = "Waived one month's of rents",
    "H2" = "Waived two months' of rents",
    "H3" = "Waived four months' of rents",
    "I"  = "Rates for residential properties",
    "I1" = "Waived up to a ceiling of $3,000 per annum",
    "I2" = "Waived up to a ceiling of $6,000 per annum",
    "I3" = "Waived up to a ceiling of $12,000 per annum",
    "J"  = "Target",
    "J1" = "Those who made public statements that call for toppling Central or HKSAR govts",
    "J2" = "Those who made public statements that criticize or call for toppling Central or HKSAR govts",
    "K"  = "Term of imprisonment",
    "K1" = "Mostly 3-5 years",
    "K2" = "Mostly 10 years to life sentence",
    "L"  = "Trial by jury",
    "L1" = "Mostly yes",
    "L2" = "Mostly no",
    "M"  = "Authority in charge",
    "M1" = "Mostly Hong Kong Police Force",
    "M2" = "Mostly National Security Agents"
)

# Specify which attribute names
attr_index <- is.na(parse_number(names(Scale_Labels)))
face_bold  <- if_else(attr_index, "bold", "plain")

# Create formulas
formula_w1 <- DF2formula(df_w1[, 4:13])
formula_w2 <- DF2formula(df_w2[, 4:17])

print(formula_w1)
print(formula_w2)

#####################
# Figure 2 (Wave 1) #
#####################

# Conjoint feature frequencies
df_w1 %>% 
    cj_freqs(formula_w1, level_order = "descending", header_fmt = "%s") %>%
    plot()

# Estimate MMs
Fig2_MM <- df_w1 %>% 
    cj(formula_w1, id = ~ID, by = ~PID, estimate = "mm", h0 = 0.5,
       level_order = "descending")

Fig2_MM_Table <- Fig2_MM %>%
    as_tibble() %>%
    select(PID, 
           Attribute = feature,
           Level     = level,
           MM        = estimate,
           SE        = std.error, 
           Lower     = lower,
           Upper     = upper, 
           p) %>%
    mutate(across(PID:Level, as.character),
           across(MM:p, ~sprintf("%.3f", .x))) %>%
    arrange(desc(PID), Attribute, Level) %>%
    mutate(Attribute = recode(Attribute, !!!as.list(Scale_Labels)),
           Level     = recode(Level, !!!as.list(Scale_Labels)))

Fig2_MM_Table

# Visualization
Fig2_MM %>%
    mutate(BY = factor(BY, levels = c("Pro-establishment",
                                      "Localist",
                                      "Pan-democrats"))) %>%
    plot(vline = 0.5, size = 2,
         xlab = "Marginal Means", header_fmt = "%s") +
    scale_y_discrete(labels = Scale_Labels) +
    guides(color = "none") +
    coord_cartesian(xlim = c(0.3, 0.7)) +
    scale_x_continuous(breaks = seq(0.3, 0.7, 0.1),
                       labels = seq(0.3, 0.7, 0.1)) +
    facet_wrap(~ BY) +
    theme(axis.text.y = element_text(face = rev(face_bold[1:36])))

#####################
# Figure 4 (Wave 2) #
#####################

# Conjoint feature frequencies
df_w2 %>% 
    cj_freqs(formula_w1, level_order = "descending", header_fmt = "%s") %>%
    plot()

# Estimate MMs
Fig4_MM <- df_w2 %>% 
    cj(formula_w2, id = ~ID, by = ~PID, estimate = "mm", h0 = 0.5,
       level_order = "descending")

Fig4_MM_Table <- Fig4_MM %>%
    as_tibble() %>%
    select(PID, 
           Attribute = feature,
           Level     = level,
           MM        = estimate,
           SE        = std.error, 
           Lower     = lower,
           Upper     = upper, 
           p) %>%
    mutate(across(PID:Level, as.character),
           across(MM:p, ~sprintf("%.3f", .x))) %>%
    arrange(desc(PID), Attribute, Level) %>%
    mutate(Attribute = recode(Attribute, !!!as.list(Scale_Labels)),
           Level     = recode(Level, !!!as.list(Scale_Labels)))

Fig4_MM_Table

# Visualization
Fig4_MM %>% 
    mutate(BY = factor(BY, levels = c("Pro-establishment",
                                      "Localist",
                                      "Pan-democrats"))) %>%
    plot(vline = 0.5, size = 2,
         xlab = "Marginal Means", header_fmt = "%s") +
    scale_y_discrete(labels = Scale_Labels) +
    guides(color = "none") +
    coord_cartesian(xlim = c(0.38, 0.6)) +
    scale_x_continuous(breaks = seq(0.3, 0.7, 0.1),
                       labels = seq(0.3, 0.7, 0.1)) +
    facet_wrap(~ BY) +
    theme(axis.text.y = element_text(face = rev(face_bold)))

##########################
# Figure 5 (Diff in MMs) #
##########################

# Combine datasets
df_combined <- bind_rows(list("Wave 1" = df_w1, "Wave 2" = df_w2), 
                         .id = "Wave") %>%
    mutate(Wave = fct_inorder(Wave))

# Preference Heterogeneity Diagnostics
## Pro-establishment
df_combined %>%
    filter(PID == "Pro-establishment") %>%
    cj_anova(formula_w1, id = ~ID, by = ~ Wave)

## Pan-democrats
df_combined %>%
    filter(PID == "Pan-democrats") %>%
    cj_anova(formula_w1, id = ~ID, by = ~ Wave)

## Localist
df_combined %>%
    filter(PID == "Localist") %>%
    cj_anova(formula_w1, id = ~ID, by = ~ Wave)

# Estimate differences in MMs (by PID)
Diff_MM1 <- df_combined %>%
    filter(PID == "Pro-establishment") %>%
    cj(formula_w1, id = ~ID, estimate = "mm_diff", by = ~ Wave, h0 = 0,
       level_order = "descending")

Diff_MM2 <- df_combined %>%
    filter(PID == "Pan-democrats") %>%
    cj(formula_w1, id = ~ID, estimate = "mm_diff", by = ~ Wave, h0 = 0,
       level_order = "descending")

Diff_MM3 <- df_combined %>%
    filter(PID == "Localist") %>%
    cj(formula_w1, id = ~ID, estimate = "mm_diff", by = ~ Wave, h0 = 0,
       level_order = "descending")

# Combining
Diff_MM <- bind_rows(list("Pro-establishment" = Diff_MM1,
                          "Pan-democrats"     = Diff_MM2,
                          "Localist"          = Diff_MM3),
                     .id = "PID") %>%
    mutate(PID = fct_inorder(PID))

Diff_MM_Table <- Diff_MM %>%
    as_tibble() %>%
    select(PID, 
           Attribute     = feature,
           Level         = level,
           `Diff in MMs` = estimate,
           SE            = std.error, 
           Lower         = lower,
           Upper         = upper, 
           p) %>%
    mutate(across(PID:Level, as.character),
           across(`Diff in MMs`:p, ~sprintf("%.3f", .x))) %>%
    arrange(desc(PID), Attribute, Level) %>%
    mutate(Attribute = recode(Attribute, !!!as.list(Scale_Labels)),
           Level     = recode(Level, !!!as.list(Scale_Labels)))

Diff_MM_Table

# Visualization
Diff_MM %>% 
    mutate(BY  = PID) %>%
    plot(xlab = "MM (Wave 2) - MM (Wave 1)", header_fmt = "%s") +
    scale_y_discrete(labels = Scale_Labels) +
    guides(color = "none") +
    facet_wrap(~BY, ncol = 3) +
    theme(axis.text.y = element_text(face = rev(face_bold[1:36])))